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ABSTRACT 

The presence of a magnetic field in the corona adds structure to the 
solar wind and almost certainly plays an important role in the energetics of 
the flow. Here I will discuss analytical and numerical modeling of gas-magne- 
tic field interactions as used to compute steady, global flow. After a brief 
and incomplete review, I will describe the approach used in, and results from 
a recent global model (Steinolfson, Suess and Wu, 1982). I will then outline 
my own ideas on the most effective ways to improve the physical content and 
numerical efficiency of these models. Throughout, I will limit myself to 
discussing solutions of the MHD equations only in order to find steady-state 
flows, even though this will often entail solving time-dependent equations. 


INTRODUCTION 

One of the more difficult problems in coronal dynamics is to self-consis- 
tently compute the large-scale interactions between the plasma and the 
magnetic field in order to find the flow geometry. The theoretical models 
that have been published generally either incorporate very important approxi- 
m-.tions in order to make the analysis tractable, or use "brute force" 
numerical solutions. As an example of an essentially analytic approach to MHD 
modeling, I describe the results from a quasi-radial flow approximation 
(Suess, Richter, Winge and Nerney, 1977; Suess, 1979; Winge and Coleman, 1974; 
Nerney and Suess, 1975; Suess, 1972). This approximation invoked flow that 
was nearly radial, deriving conditions for this requirement as a result from 
the analysis. Single-fluid, polytropic flow with no dissipation was assumed, 
and the flow was taken to be axisymmetric. The geometrical approximation 
means that the analysis cannot be applied near a magnetic cusp. The model was 
applied to the northern polar coronal hole of 1973 because observations 
provided a density distribution throughout the hole and showed that the hole 
was essentially axisymmetric (Munro and Jackson, 1977). The observed density 
was matched throughout the entire volume, and the observed geometry of the 
boundary was matched to a streamline, minimizing the impact of the polytrope 
assumption. The results are shown in Figures 1 and 2. Important deductions 
include an "effective" temperature of 1-2 x 10® deg., a field strength of 
0.5 to 1.0 gauss, and a flow speed of a few up to 150 km/s at 2 solar radii. 
By 5 solar radii, the flow speed varies from less than 50 to over 300 km/s 
while the magnetic field is approximately uniform across the hole. 
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Figure 1. The observed boundary of the northern polar coronal hole (solid curve, 
from Munro and Jackson, 1977) and streamlines beginning at 2 solar radii and 
polar angles of 20 and 50 degrees. The streamlines have been computed using a 
quasi-radial approximation to the MHD equations of motion together with latitude 
dependent boundary conditions (from Suess, , 1977). 
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Figure 2. Latitudinal variations of the radial magnetic field, temperature and 
radial velocity at 2 and 5 solar radii. The radial magnetic field, temperature 
and observed density variation at 2 solar radii are the boundary conditions for 
this, and the streamline results in Figure 1 (Suess, _et £l . , 1977). 
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Significant energy deposition is implied between 2 and 5 solar radii. 
However, because this analysis cannot treat global flows that include regions 
of closed magnetic fields, it cannot be carried much beyond this particular 
example. The only obvious extensions would be to relax the polytropa 
assumption to permit thermal dissipation, and to include momentum and energy 
source terms. 

In a corona with both closed, magnetostatic regions, and open, coronal- 
hole-like regions, the problem of the global gas-magnetic field interactions 
is both totally nonlinear and geometrically complex, almost requiring numeri- 
cal solution. Probably the most pioneering work was that by Pneuraan and Kopp 
(1971), who constructed a solution for an axisymmetric , isothermal streamer 
configuration on a two-dimensional grid. Their solution was found through a 
procedure which iterated on the current distribution, and they assumed 
isothermal flow. 

The ultimate fieldline geometry computed by Pneuman and Kopp is shown in 
Figure 3. Originally, there were questions regarding their assumption about 
the nature of the cusp. However, these were laid to rest using a numerical 
solution of the time-dependent equations of motion (Endler, 1971; Weber, 
1978). If the situation is treated as an initial-boundary value problem, then 
it is guaranteed that the final state is a true equilibrium configuration. In 
this approach, one starts with an initial state consisting of an essentially 
arbitrary choice for the fluid and magnetic field variables. The numerical 
solution of the time-dependent equations, starting with this non-equilibrium 
state, then gradually approaches a stationary coronal configuration. In 
applying this approach, it was not only shown that the cusp geometry assumed 
by Pneuman and Kopp was a valid equilibrium configuration, but also that the 
solution was stable. 

These results demonstrated that treatment of steady, global coronal flow 
as an initial-boundary value problem is: (i) feasible, (ii) efficient, and 
(iii) probably the most powerful approach to the general problem. I would now 
like to describe results from an application of this approach to polytropic 
rather than isothermal flow, for a survey of the effect of varying magnetic 
field strength (Steinolfson, et. al, 1982). 


THE PROBLEM AND THE EQUATIONS 


The assumptions are that the corona can be described by polytropic, 
axisymmetric, single fluid flow. The initial state consists of a hydrodynamic 
solution to the steady radial equations of motion for a polytropic gas, 
superimposed on a dipole magnetic field. No explicit dissipation is included. 
This is a generalization over previous work by relaxing the assumption of 
isothermal flow. With these assumptions, the time-dependent equations of 
motion, in MKS units, are: 
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where the dependent variables are the density, P, radial velocity, Vr. 
meridional velocity, v^, pressure, p, radial magnetic field, Br and meridion- 
al magnetic field, The Independent variables are the radius, r, and the 

colatitude, 0. The constants are the polytropic index, y» magnetic permeabil- 
ity, p, universal graviational constant, G, and the solar mass, M^. 

These equations are solved between 1.0 and 5.0 solar radii, and from the 
pole to the equator - the solution is symmetric about the equator. The grid 
spacing is 0.1 solar radii in the radial direction, 2.5 degrees in the 
meridional direction, the numerical solution uses a modified Lax-Wendroff 
^iff*®i“®ncing scheme, and the time step is chosen to be the maximum allowable 
from the usual stability criterion for dissipationless schemes; i.e.. 


At = min (At , At ) 
r 0 


( 2 ) 


where 

At = r/A 


r r 

and 

Atg = r0/X^ 


and maximum eigenvalues (the sum of the fluid velocity 

and the characteristic velocity) in the radial and meridional directions. A 
smoothing term is used to reduce numerical oscillations, and it was necessary 
to watch for nonzero values of div.B, which otherwise often exist in numerical 
solutions of multi-dimensional time-dependent MHD problems (Brackbill and 
Barnes, 1980). At the inner boundary, two of the six radial characteristic 
directions are negative, and consequently, information from the region of 
interest propagates upstream to the boundary. In this case, four dependent 
variables at the lower boundary can be specified arbitrarily, and two must be 
calculated from some form of compatibility relations (Steinolfson and Nakaga- 
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length is proportional to the 
velocity magnitude (Steinolf- 
son, et al. , 1982) • 







wa, 1976). Strictly speaking, the compatibility relations are equations that 
can be derived from the equations of motion which must be satisfied by the 
dependent variables in each of the characteristic directions. Steinolfson and 
Nakagawa have shown that first-order or second-order (linear) extrapolation 
often works as well as using the more canplex compatibility relations. 
Steinolfson, et . al used linear extrapolation to obtain the pressure and 
meridional magnetic field at the inner boundary. The radial magnetic field 
was held at its initial value so that the total magnetic flux through the 
solar surface remains constant. The radial velocity was also held constant 
with the exception that it could decrease to zero at the surface inside the 
closed field region. The meridional velocity was calculated so that the total 
velocity and magnetic field were parallel at the surface. The surface density 
was selected so that p/p’'^ was constant. 

In the initial state of spherically symmetric flow and a dipole potential 
field, the boundary values at 1 solar radius, and the system constants are: 

= 1.8 X 10^ degrees 

n = 2.25 X 10^ cm ^ 
o 

Y = 1.05 


B = 


2n k T 
o o 

o 


(3) 


B = |b| at 1 solar radius, at the equator. 

% 

Using this definition of the plasma beta, results are summarized below for 
values of beta from 0.1 to 100. Values for the reference magnetic field of 
0.83 G and 2.35 G at the equator yield values for the plasma beta of 4 and 
0.5, respectively. 


RESULTS: RELAXATION TO A STEADY STATE 

Beginning with the previously stated initial configuration of a dipole 
potential field superimposed on a spherically symmetric flow, the relaxation 
then precedes in time until the solution is approximately steady - meaning the 
solution did not change appreciably over a period of 2 hours. The interme- 
diate states themselves have little physical meaning. However, the relaxation 
time to a steady state does have physical meaning - it is typical of the time 
the corona would take, given the assumptions of the model, to return to 
equilibrium after a large-amplitude perturbation. Figure 4 shows the evolu- 
tion of the coronal magnetic field (solid lines) and velocity (vectors) for 
beta = 0.5. The initial state is in panel (a), and subseqent states at 4, 8 
and 12 hours are shown in panels (b), (c) and (d) respectively. There is 
little change after 12 hours. It is easy to see the field lines evolving from 
a closed dipole field to a coronal streamer with the closed field lines lying 
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beneath and adjacent to open jld lines. The dashed line is the sonic curve 
and the dotted line is the Alfven curve. The sonic curve is displaced inwards 
in the final state, except for a small region around the equator, due to the 
general increase in the velocity. In the closed region, the pressure and 
density are increased over their initial values, and the velocity is 
approximately zero. Figure 5 shows cuts through the final configuration, at 
the pole and the equator. The high density and pressure in the closed region 
are evident, as is the zero velocity. Enhanced flow in the center of the open 
region is due essentially to the effects described by Kopp and Holzer (1976). 
However, the overall results are not exactly equivalent because they held the 
energy per gram, or flow speed infinitely far from the sun fixed, whereas here 
the boundary condition on temperature and velocity was held fixed on the sun. 
The Kopp and Holzer study dealt with the effects of changing "spreading 
factor" (overall flux tube divergence), whereas here a change in divergence 
occurs naturally along different streamlines. 

Some of the effects of differing divergence can be simulated by changing 
the magnetic field strength. This was done by changing beta in steps from 
100.0 to 0.1. Again, the results cannot be directly compared with those of 
Kopp and Holzer because the energy per gram at infinity is not a constant. 
The results from changing field strength are summarized in Figures 6 and 7. 
Figure 6, showing the maxiraim pressure and density in the closed region after 
scaling by the initial flow values, indicates the physical relevance of the 
beta parameter. The curves remain near a value of unity until beta < 1.0. 
Then, they rise steeply for values of beta less than ca. 0.5. At the same 
time, there is a corresponding general increase in the height of the streamer. 
However, the width of the base of the streamer has only a weak dependence on 
beta. 


Figure 7 shows the latitudinal variation of the pressure, density and 
temperature at five solar radii, for a range of beta on either side of unity. 
Again, these variables are scaled by their initial values. The velocity is 
shown in Figure 8. For values of beta >= 1 , the only substantial effect is 
that the velocity is slightly depressed and the density is enhanced in the 
streamer. As beta is reduced below unity, a more rapid change begins to 
appear over the open region. Simultaneously, the velocity in the hole begins 
to increase rapidly, and the density and temperature begin to decrease. This 
is, in general terms, the same phenomenon observed by Kopp and Holzer (1976) 
and Steinolfson and Tandberg-Hanssen (1977), and is due to the divergence 
along a streamline being more than the underlying radius-squared divergence - 
which is commonly called a spreading factor of more than unity. 

Note that the maximum velocity at 5 solar radii is not at the center of 
the hole, but instead near the edge at a polar angle of 60 degrees at 5 solar 
radii. This reflects a combination of the boundary condition on velocity 
being spherically symmetric, together with the overall spreading factor 
between 1 and 5 solar radii being dominated by potential field-like effects 
near the sun. The spreading factor is largest near the edge of the hole, 
although it is also larger than unity near the center of the hole. As beta 
decreases below unity, the increasing spreading factor causes the sonic 
critical point to begin moving rapidly inward (from 3.5 to 2.9 solar radii for 
beta decreasing from 0.5 to 0.1). As in earlier studies, the consequence is a 
rapid increase in local flow speed (but not the terminal flow speed), with the 
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Figure 5. Radial distribution of the density, pressure and radial velocity ini- 
tially (dashed line) and in the coronal streamer for beta = 0.5, at the pole and 
equator (Steinolf son, je;t al . , 1982). 



BETA 

Figure 6. The dependence of the maximum pressure and density in the magnetically 
closed region on the reference beta (Steinolf son, ^ al. , 1982). 
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largest increase occurring where the spreading factor is largest. 

The result of this simulation and survey of magnetic field stength 
effects has been to produce a self-consistent streamer-coronal hole configura- 
tion. The streamer, in which the closed field region extends slightly more 
than 1 solar radius above the surface and whose width is somewhat less than 
the height, is quite similar to a typical streamer in its general features. 
The coronal hole is geometricaly similar to, for example, the polar coronal 
hole observed during the Skylab period , and the density at the center of the 
hole is less than at the edges. However, there is good reason to believe the 
velocity should also increase towards the center of the hole (Suess, et. al , 
1977). That it doesn't is due to the choice of boundary conditions and the 
lack of extended energy and momentum addition other than what is implicit in 
the polytrope formulation. The simulation shown in Figures 1 and 2 gave a 
detailed fit to the polar coronal hole data reported by Munro and Jackson 
(1977). There, it was shown that a large variation in density and temperature 
was required at the inner boundary - 2 solar radii - in order to fit the 
observations. In order to fit those data using the model of Steinolfson, et 
al . , it would be necessary to invoke a corresponding variation at 1 solar 
radius, or to deposit energy into the flow between 1 and 2 solar radii. The 
conclusion then must be that to simulate at least one specific coronal hole, 
energy and momentum must be deposited in the flow between 1 and 2 solar radii 
because there is little evidence for the alternative possibility of signifi- 
cantly higher temperatures and densities beneath the hole, at the transition 
region level. 

Before leaving this topic, there are two additional important specific 
results worth noting (R. S. Steinolfson, pers. comm.). First, there has 
been a question of whether the final configuration shown in Figure 4 is 
unique. To address this first question, Steinolfson began with an initial 
state that was radically different that that shown in Figure 4. It had the 
same spherically symmetric flow field, but now with a strictly radial magnetic 
field whose strength at the surface varied in the same way as the radial 
component of a dipole field. This exercise requires that field lines can 
"reconnect" as they are advected inward through the outer boundary in order to 
form the closed field region of the streamer. This is possible here not 
because ohmic diffusion is specifically included, but because there is 
sufficient nwnerical diffusion to allow a similar process to occur. The 
simulation produced exactly the same final configuration as was found when 
starting with a dipole potential field, so that the conclusion is: the present 
example has a unique solution. The second question has to do with the 
suggestion that standing shock waves may exist for certain ranges of spreading 
factor and initial conditions (Holzer, 1977; Habbal and Tsinganos, 1983). The 
study of Steinolfson, et cannot completely answer the question of whether 
these shocks ever actually occur in flows where the transverse pressure 
balance in the presence of a magnetic field is explicitly computed, because a 
beta of 0.1 is not sufficiently small to produce the extreme spreading factors 
invoked in some of those studies. However, these computations for beta 
between 100 and 0.1 are probably sufficient to cover most cases of relevance 
to the sun. Since no standing shocks were observed, it is suggested that the 
phenomenon may not be of importance in the solar atmosphere. 
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PROSPECT FOR FUTURE MODELING 


The most interesting direction that can be taken in the development of 
numerical models of the steady corona would be to include a more detailed 
description of the energetics of the flow. The polytropic description is 
limited by many well known deficiencies. An obvious step would be to treat 
thermal conduction explicitly, and later to expand the calculation to deal 
with non-collisional forms of the conductive energy transfer. Other similar, 
but less important generalizations would be to include radiative losses and 
ohmic diffusion. For detailed models of specific observations, until the 
actual processes are discovered, it may also be necessary to include an 
empirical description of energy and momentum deposition. As stated earlier, 
the model described by Steinolfson, et al. used a modified Lax-Wendroff 
differencing scheme, which incorporates explicit time-differencing. Because 
of the Courant condition for stability on the maximum time-step, such models 
are effectively unuseable for calculations including thermal or ohmic diffu- 
sion because the maximum allowable time step decreases linearly with the 
largest characteristic speed. Literally thousands of time-steps would be 
required with thermal diffusion in order to do the analogous problem to the 
relaxation described earlier. 

The alternative is to use implicit time-differencing vdiich, although 
algebraically complex, is not limited by the Courant condition or round-off 
errors proportional to (At/Ax) ^ - as is the case with the Dufort-Frankel 
method (Richttneyer & Morton, 1967). In fact, it is probable that these 
schemes would result in a considerable improvement in computational efficiency 
- even for the problem treated by Steinolfson, et al . The reason for this is 
that, because no particular physical interest exists for the intermediate 
steps between the arbitrary initial state and the final steady-state, time 
steps too large to resolve specific intermediate dynamic fluctuations can be 
used and still arrive at the desired final state. This is because the 
steady-state solution, if one exists, is found totally independent of the step 
size. The one qualification is that if multiple steady-state solutions exist, 
then there is the possibility that the solution that is actually realized will 
depend on the step size because a sufficiently large dynamic fluctuation can 
cause the asymptotic approach to a final solution to jump to an alternative 
branch. 

To the beat of my knowledge, no implicit scheme has been applied to the 
problem described by Steinolfson, et As a substitute, I will summarize 
the results from a prototype calculation on one-dimensional, transonic, 
thermally conductive flow (Suess, 1982) , and outline how the analysis is 
extended to multi-dimensional magnetohydrodynamic flow. First, for the 
one-dimensional problem, the equations are: 


dt 
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Where p is the density, is the radial velocity, p is the pressure, e Is 
the total (kinetic plus internal) energy per unit volume, and M is the mass 
of the Sun. 

p = 2nkT=2p{ J I r s 2pRT, 
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The details of the calculation are described by Suess (1982). The solution 
uses three-level implicit time-differencing that is unconditionally stable for 
any time-step. Solutions are found on a variable temporal and spatial grid 
using a near-conservation form of the equations and occasionally with a small 
amount of artificial dissipation inserted in such a way as to preserve the 
formal accuracy of the differencing scheme. The spatial grid is not split, 
and the solution is fully implicit - treating all equations simultaneously. 
The algorithm uses a Taylor series expansion and Jacobian matrices to deal 
with nonlinear coefficients so that there is no iteration of any kind between 
time levels. Finally, although three levels in time are used, the algorithm 
requires only two levels to be stored - thereby minimizing memory require- 
ments. The algorithm is applied here to a one-dimensional problem, but I 
will illustrate, in general terms, application to multi-dimensional problems 
after the example. The example is a demonstration of the efficiency of such 
algorithms on a classical problem. The initial state is isothermal solar wind 
flow with_a temperature and density at 1 solar radius of 2x10^ degrees and 
lO"*"' cm ^ respectively, resulting in a flow speed at 1 AU of £a. 400 km/s. 

Figure 9 shows the relaxation at several different radii, beginning with the 
third time-step. The time step is allowed to grow if the solution is changing 
sufficiently slowly. The relaxation precedes rapidly at first, with some 
large oscillations and overshoot, settling down to a steady state inside 30 
solar radii after 75-100 hours. The final time step was 22.2 hours and 19 
steps were required for a total elapsed time of 108 hours. 

The relaxation provides the correct, known solution and the relaxation 
time is essentially the advection time. This example is approximately equal 
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in computational efficiency to solving the steady-state equations in the 
standard way. As soon as an additional process such as radiative loss is 
introduced, the relaxation becomes the most efficient technique. 


To indicate the procedure for multi-dimensional flow, I begin with the 
form of the^time-d|fferenced equations in the one-dimensional case: 

(5) 
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where spatial differences have not yet been introduced, L is a matrix 
operator, AU” is the vector of variables in conservation form - being the 
unknown at the next time step, and F is a matrix of known quantities at the 
present and previous time steps. Solution involves a matrix inversion or 

solution of a system of equations in order to find AlP^ . In multiple 
dimensions, an approximate factorization can be made to give a corresponding 
equation of the form (for two dimensions): 

L- |~> = F + cross terms. (6) 

L 3r J 80 


which is solved in two steps: 

I- 2 T 
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cross 


AU" = AU^* 


terms. 


That is, intermediate solutions are found by sweeping in one direction through 
the grid, at eac’.i time step, and the final solution for that time step is the 
produced by sweeping in the other direction through the grid. This is 
possible because the equation has been "factored" into operators depending on 
spatial gradients in only one direction. In the presence of magnetic fields, 
the factorization is not complete - some cross terms involving only the 
magnetic field remain. These cross terras are treated explicitly. 

This factored implicit algorithm has been applied to a variety of 
ordinary hydrodynamic problems (Beam and Warming, 1976), with great success. 
It has been found generally to have superior stability limits and great 
flexibility. There is every reason to believe that it would work equally as 
well for determining steady-state structure in the corona. 


CONCLUSION 

I have described the approach to global modeling of the corona as an 
initial-boundary value problem in which the configuration is allowed to relax 
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in time from some initial state into a steady state. It is established that 
the technique is powerful, flexible and stable. The results include a survey 
of changing magnetic field strength effects and address questions about the 
uniqueness of the steady state and the possibility of standing shock waves. 
Application of advanced numerical methods now holds promise for being able to 
extend the present results by addressing the all-important problem of energy 
transfer in the corona, even in the presence of magnetic fields. 
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